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Summary 

The Integrated Force Method has been developed in recent 
years for the analysis of structural mechanics problems. This 
method treats all independent internal forces as unknown 
variables that can be calculated by simultaneously imposing 
equations of equilibrium and compatibility conditions. In this 
paper a finite element library for analyzing two-dimensional 
problems by the Integrated Force Method is presented. 
Triangular- and quadrilateral-shaped elements capable of mod- 
eling arbitrary domain configurations are presented. The ele- 
ment equilibrium and flexibility matrices are derived by 
discretizing the expressions for potential and complementary 
energies, respectively. The displacement and stress fields within 
the finite elements are independently approximated. The dis- 
placement field is interpolated as it is in the standard displace- 
ment method, and the stress field is approximated by using 
complete polynomials of the correct order. A procedure that 
uses the definitions of stress components in terms of an Airy 
stress function is developed to derive the stress interpolation 
polynomials. Such derived stress fields identically satisfy the 
equations of equilibrium. Moreover, the resulting element 
matrices are insensitive to the orientation of local coordinate 
systems. A method is devised to calculate the number of rigid 
body modes, and the present elements are shown to be free of 
spurious zero-energy modes. A number of example problems 
are solved by using the present library, and the results are 
compared with corresponding analytical solutions and with 
results from the standard displacement finite element method. 
The Integrated Force Method not only gives results that agree 
well with analytical and displacement method results but also 
outperforms the displacement method in stress calculations. 


Introduction 

The finite element stiffness method, which is based on an 
assumed displacement field, has become the method of choice 
for solving a wide variety of problems in structural mechanics. 
The advantages of the stiffness method include (1) the capabil- 
ity to efficiently and accurately model domains with complex 
geometric configurations and varying material properties and 
(2) the capability to accurately analyze problems with geo- 
metrical and material nonlinearities. The development of finite 
stiffness elements and their corresponding formulations has 
been a subject of extensive research, much of which has been 
summarized in textbooks such as references 1 to 3. 

Shortcomings of the assumed displacement method have 
been observed in the analyses of certain classes of problems, 
such as modeling nearly incompressible materials, bending of 
thin plates, and optimizing structures (refs. 4 and 5). Moreover, 
since stresses are calculated indirectly by using displacement 
derivatives, the accuracy of stress predictions may be reduced. 
Two alternative finite element formulations may be utilized to 
analyze the aforementioned problems and to calculate stress 
more accurately: (1) the hybrid stress method (refs. 6 to 8), and 
(2) the force method (refs. 9 to 11). Because both of these 
formulations have certain disadvantages compared to the 
assumed displacement method, their use and availability in 
general purpose programs has been limited. In the hybrid 
method, the flexibility matrix must be inverted in order to 
generate the element stiffness matrix; this can become a com- 
putational burden, especially if high order approximations of 
stress fields are required. In the standard force method, on the 
other hand, an auxiliary statically determinate structure and a 
corresponding set of redundant forces must be selected. This 



process is not easily adapted to computer automation. Several 
attempts have been made to improve the process by which 
redundancies are selected. The pertinent formulations were 
summarized by Kaneko et al. (ref. 5). All of these procedures, 
however, either resulted in matrices with certain undesired 
properties or lacked a physical interpretation, which made them 
unattractive to the engineering community and led to the 
demise of the standard force method. 

An alternate formulation, termed the IntegratedForce Method, 
has been developed in recent years to analyze problems in 
structural mechanics (refs. 12 to 15). In the Integrated Force 
Method all independent forces, not just the redundants, are 
treated as unknown quantities that can be calculated by simul- 
taneously imposing both equilibrium and compatibility condi- 
tions. Procedures have been developed (refs. 16 to 19) for 
generating compatibility conditions that yield sparse and banded 
matrices and can be easily adapted to computer automation. 
The initial applications of the Integrated Force Method to static 
analysis (ref. 20), vibration analysis (ref, 21), and optimization 
of trusses (ref. 22) have shown that the Integrated Force Method 
has certain advantages over the displacement method, both in 
accuracy and computer efficiency. 

This study presents formulations to develop finite elements 
for two-dimensional structural analysis and a comprehen- 
sive finite element library of two-dimensional elements. Both 
triangular- and quadrilateral-shaped elements capable of mod- 
eling arbitrary configurations of the domains being analyzed 
are considered. The displacement and stress fields within an 
element are independently approximated. The displacement 
field is interpolated by using the functions employed in the 
standard displacement method. Stress fields are approximated 
by using complete polynomials of the appropriate order, whose 
coefficients are unknown independent forces. The equations 
describing the components of the stress tensor can be derived 
from the Airy stress function for an element, which is written 
in terms of a complete polynomial of a certain order. The 
resulting stress fields identically satisfy the equations of equi- 
librium. The element matrices generated with these stress fields 
are not sensitive to the orientation of the element’s local 
coordinate system. A method to calculate the number of zero- 
energy modes is also developed, and the present elements are 
shown to be free of spurious zero-energy modes. The effect of 
reducing the number of the element’ s independent forces is also 
investigated. 

To establish the validity of the elements and to assess then- 
relative performances and compare the Integrated Force Method 
with the well established displacement method, the present 
library is used to solve a variety of problems in two-dimensional 
elasticity. The results obtained with these elements are also 
compared with the corresponding analytical solutions, and 
there is good agreement. For stress calculations, the Integrated 
Force Method performs better than the standard displacement 
method. 


Development of the Finite Elements 

Governing Equations of the Integrated Force Method 

The governing equations of the Integrated Force Method are 
briefly presented here in order to introduce the notation. Sym- 
bols used are defined in appendix A. (A detailed description of 
the formulation can be found in refs. 12 to 14.) Finite elements 
are used to discretize a continuous object, which then has N t 
displacement degrees of freedom and m independent forces. In 
the Integrated Force Method all independent forces represent 
unknown quantities, not just the redundants, as is the case in the 
standard force method (refs. 9 to 1 1). The unknown forces are 
obtained from the following sets of equations: 

[B m ){F} = {P} (1) 

which represents n equations of equilibrium, written for nodes 
where displacements are unknown, and 

[C][G]{F} = {8 0 } (2) 

which represents r~m-n compatibility conditions. 

Here n = N t ~N s ; N s is the number of prescribed displacement 
degrees of freedom; {F} is the m-c omponent vector of un- 
known independent forces; [BJ is the n x m part of the system 
equilibrium matrix corresponding to the nodes where external 
loads are prescribed; [G] is the m x m system flexibility matrix; 
[C] is the rxm compatibility matrix; {P} is the n-component 
vector of equivalent nodal loads; and {5 0 } is the r-component 
effective deformation vector, which is calculated as 

{5,} = -rC]{p 0 } (3) 

where { %} is the vector of initial deformations. The generation 
of compatibility conditions is described in detail in references 
16to 19. Sets of equations ( 1 ) and (2) can be combined to obtain 
the system of equations for unknown forces as 


where 


[S]{F} = {P*} 


( 4 ) 


[S] = 


[C][G] 


and 



( 5 ) 


After the force vector is calculated from equation (4), the vector 
of unknown nodal displacements {U} can be obtained as 

{U} = [J][G]{F} (6) 
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where [ J] is the n x m deformation matrix that represents the top 
n rows of the transpose of the matrix [S]" 1 . The vector of 
unknown support reactions {R} can be calculated from 

{R} = [BJ = {F} (7) 

where [BJ is the portion of the system equilibrium matrix that 
corresponds to nodes with prescribed displacement boundary 
conditions. 

Element Matrices 

From equations (1) and (2) we see that the system of 
equations for the unknown forces consists of two sets of 
relations: (1) equations of equilibrium and (2) compatibility 
conditions, which can be expressed in terms of forces by using 
the strain-stress law. These sets of relations are first established 
on the element level, and then the assembly procedure (refs. 13 
and 14) is used to derive the system given in equation (4). 

Equilibrium equations and deformation-force relations for a 
finite element may be written as 

{P e } = (BJ{F e } (8a) 

® e } = [G e ]{F e } (8b) 

where {P e } is the vector of equivalent nodal forces for the 
element e\ {F^} is the element vector of independent forces; 
[B e ] and [G e ] are element equilibrium and flexibility matrices, 
respectively; and {$ e } is the vector of element deformations. 
Note that equation (8b) represents the discretized constitutive 
relations for the element. The components of the vector {f3 £ } 
are the generalized deformations that correspond to internal 
forces {F^}. 

The expressions for the two element matrices can be derived 
by using the expressions for potential and complementary 
energy, respectively. In the Integrated Force Method, indepen- 
dent displacement and stress interpolations are employed to 
give 

{u} = [N]{U e } (9) 

{a} = [Y]{F e } (10) 

Here {u} T = {uv} and { o } T ={G x CFyTxy} are the displacement 
and stress vectors at a location within the element; {U e } is the 
vector of element nodal displacements; [N] is the matrix of 
displacement interpolation functions; and [Y] is the stress 
interpolation matrix. The strain vector, {e} T = {e x Sy Yxy L is 
obtained by differentiation of the displacement field and is 
given as 


{£} = [Z]{iy (11) 

where fZ] = [L] [N], and [L] is the matrix of differential 
operators. 

Equilibrium matrix . — The expression for the element equi- 
librium matrix can be obtained from the strain energy A p of the 
element: 

A p =|j v {e} T {0}rfV (12) 

where V denotes the domain of the element in the discrete form. 
Substituting equations (10) and (11) into equation (12) yields 
the strain energy A p expressed as 

A p =|{U e } T [B e ]{F e } (13) 

where the equilibrium matrix [B e ] is 

[B g ] = J v [Z] J [Y]dV (14) 

Flexibility matrix . — -The expression for the element flex- 
ibility matrix can be obtained from the complementary energy 
A c of the element: 

A c=\\ v M T [DKo}dV (15) 

where [D] is the compliance matrix of the material, in the 
discrete form. Substituting equation (10) into equation (15) 
gives A c as 

A c = ±{F e } J [G e W e } (16) 

where the element flexibility matrix [G e ] is 

[G e ] = J y [Y] t [D][YW (17) 

Stress Field Approximations 

The approximation of the stress components and the con- 
struction of the stress interpolation matrix [Y] is discussed in 
this section. Note that inequations (14) and (17) thematrix [Y] 
appears in both the equilibrium and the flexibility matrix 
definitions. It is, therefore, important to properly devise stress 
interpolation polynomials in order to obtain accurate results. In 
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this study, a method was developed that uses an Airy stress 
function given in terms of complete polynomials. 

The Airy stress function 0 for a location (x, y ) within an 
element can be written as a complete polynomial of order p: 

P 

0(x,y) = ^ CjX p ~ j y j (18) 

7—0 

where Cj, for j = 0,1,2,...,/?, are constants, and x and y are 
Cartesian coordinates of the point in the local coordinate system 
of the element. The local coordinate systems for various ele- 
ment shapes are depicted in figure 1. The components of the 
stress tensor can be obtained by using the definition of the stress 
function (ref. 23) as follows: 


Rewriting equations (19) yields 
P-2 

^=X^~ V (20a) 

7=0 

P-2 

(20b) 

7=0 

P-2 

(20c) 

7=0 


P-2 

^ = X c ' +2 ° + 1)0 + 2)xP ~ 2 ~ iyj (i9a) 

7=0 

P-2 

CF y = - j)(p~j - l)x p - 2 ~ j y j (19b) 

7=0 

P-2 

V=-Xrf C i+lO‘ + l)(P-J-l)^ P-2 ~V (!9c) 

7=0 
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Figure 1 . — Finite element library for Integrated Force Method. 


Now the coefficients of the polynomials in equations (20) can 
be considered element forces F for i = 1,2,. ..,3 (p - 1), and F. 
can be expressed in terms of (p + 1) constants Cj : 

F.~<k(c 0 ,C 1 ,...,C p ) for i = l,2,,..,3(p — l) (21) 

where 0/ are linear functions of constants Cj. Thus, not all forces 
F are linearly independent. Final stress field interpolation 
polynomials can be obtained by eliminating the dependent 
forces; this results in (p + 1) independent forces when the stress 
function 0 is written as a complete polynomial of order p. Such 
stress fields are complete polynomials of order (p- 2). Expres- 
sions for the stress fields used in this study are provided in 
appendix B. 

The preceding procedure can be demonstrated by deriving 
the linear terms. For this case, the stress function is represented 
as a complete cubic polynomial: 

0{x, y) CqX^ + C^x^ y + C 2 xy^ + C^y ^ ( 22 ) 

Substituting p = 3 into equations (19) yields the following 
expressions for the stress components: 


<T® = 2C 2 y + 6C 3 y = F^x 4- F 2 y 

(23a) 

<7® = 6C 0 x + 2C±y ■■= F 3 jc + F^y 

(23b) 

t^ = -2C 1 x-2C 2 y = F 5 x+F 6 y 

(23c) 


where the superscript 1 denotes the linear terms in the stress 
polynomials. From equations (23) we see that six coeffi- 
cients F t are expressed in terms of four independent constants 
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Cj\ thus, only four E forces are linearly independent. By 
eliminating two forces from equations (23), we have obtained 
the linear stress terms given in equations (B2). The constant, 
quadratic, and cubic terms can be derived similarly, and the 
corresponding stress terms are given in equations (Bl), (B3), 
and (B4), respectively. Stress field representations in terms of 
a complete polynomial of order p can be obtained by combining 
the expressions for orders 0,1,2,...,/?. The resulting cubic stress 
field interpolation, which is given in equations (B5), contains 
1 8 independent forces. The constant stress field can be obtained 
by retaining the first 3 forces in equations (B5); for the linear 
interpolation, the first 7 forces are retained, and for the qua- 
dratic interpolation, the first 12 forces are retained. The stress 
fields given in equations (Bl) to (B4) identically satisfy the 
equations of equilibrium at any point inside the domain of the 
element. The resulting element matrices have the correct rank 
for arbitrary orientation of the element’s local coordinate axes. 
They are also invariant with respect to coordinate 
transformation (ref. 4). 

Approximations with complete polynomials may yield indi- 
vidual elements with a large number of independent forces. 
Moreover, as Spilker and Singh (ref. 24) observed, in hybrid 
method applications, high order stress field approximations 
may lead to overly rigid models. Thus, it may be necessary to 
reduce the number of independent forces in stress field repre- 
sentations while preserving all the desired properties of the 
resulting element matrices. The compatibility condition 

V 2 «r x + <J y )=0 (24) 

(suggested in ref. 24) can be imposed to reduce the number of 
independent forces. Note that equation (24) is identically 
satisfied for stress fields represented by zero- and first-order 
polynomials. By applying equation (24) to quadratic and cubic 
terms, we obtain reduced quadratic and cubic polynomials, as 
given in equations (B6) and (B7), respectively. 

Spurious Zero-Energy Modes 

The stress fields given in equations (B 1) to (B7) were derived 
without any reference to the shape or the number of kinematic 
degrees of freedom of a considered element. The number of 
independentforces, however, couldnot be chosen arbitrarily. The 
number of kinematic degrees of freedom n e and the number of 
independent forces m e for element e must satisfy the relation 
m e >n e - 1 (refs. 25 and 26), where / is the number of rigid body 
modes of the element. Pian and Chen (ref. 26) showed, however, 
that in the application of the hybrid method this condition is only 
necessary, not sufficient, for the element matrices to have the 
correct rank. They also devised a technique, based on energy 
considerations, to detect spurious zero-energy modes, and they 
developed a means to suppress them. Spilker et al. (ref. 4) have 
shown that approximating stresses with complete polynomials of 


the appropriate order produces the correct rank. The methodology 
of Pian andChenis dsedherein to show that in the IntegratedForce 
Method the correct rank of the element equilibrium matrix 
ensures the absence of spurious zero-energy modes. 

The expression for internal energy A c can be rewritten by 
substituting [D]{a} = {e} into equation (15) to obtain 

A c =±j v {o} r {e}dV (25) 

Substituting equations (10) and (11) into equation (25) gives 
the internal energy written as 

A C =~{F} T [B] T {U} (26) 

From equation (26) we can see that if the element equilibrium 
matrix has the correct rank, there are only / zero-energy modes 
present that correspond to rigid body modes of the element. 
Thus, spurious zero-energy modes can be eliminated by con- 
structing stress fields such that the resulting equilibrium matrix 
has the rank n r >n e - 1. 

Element Library 

Stress fields derived in previous sections can now be used to 
develop a comprehensive finite element library for two- 
dimensional stress analysis by the Integrated Force Method. 
Let us consider both triangular and quadrilateral elements. 
Isoparametric functions (ref. 27) can be employed in equa- 
tion (9) for both types of elements. For stress field approxima- 
tions, the element’s local coordinate systems Oxy can be 
defined such that the origin coincides with the centroid of the 
element, and the local coordinate axes are parallel to the global 
axes. Such an orientation avoids rotation of the coordinate sys- 
tems, saves CPU time, and does not affect the response when 
stress fields with complete polynomial approximation are used. 

The element library is depicted in figure 1. It includes two 
elements developed by Nagabhushanam (J. Nagabhushanam, 
Indian Institute of Science, Bangalore, India, personal commu- 
nication, 1992) and an element suggested by Pian (ref. 6), 
which has four nodes and five independent forces. These 
elements are also implemented for comparison purposes. The 
element names employed here consist of three parts: the first 
three characters describe the shape of the element, the next two 
digits denote the number of element nodes, and finally, the 
number following the underscore indicates the number of 
independent forces used in the interpolation of the stress field. 
Features of the present elements are enumerated in the follow- 
ing sections. 

Three-node triangular elements: TRI03_03, TRI03_05, 
and TRI03JQ7. — Three-node triangles have six displacement 
degrees of freedom; thus three independent forces are necessary 
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to ensure the correct rank of element matrices. Three different 
stress fields were implemented. The constant stress field was 
used for element TRI03_03 . A complete linear stress field was 
implemented for element TRI03_07. And the stress field given 
in equations (B8), which is used extensively in the hybrid 
method (refs. 6 and 11), was used for element TRI03_05. Note 
that element TRI03_03 contains the minimum number of 
independent forces. Such elements will be referred to as stati- 
cally determinate elements. 

Six-node triangular elements: TRI06_09 > TRI06J11, and 
TRI06_12— Six-node elements have 12 displacement degrees 
of freedom; thus 9 independent forces are necessary in the 
stress field approximation. This allows quadratic polynomials 
to be used for stress interpolations. Complete quadratic polyno- 
mials with 12 independent forces were used for element 
TRI06_12, and reduced quadratic polynomials with 11 inde- 
pendent forces were used for element TRI06_1 1 . For compari- 
son the stress field suggested by Nagabhushanam (Indian 
Institute of Science, Bangalore, India, personal communication, 
1992; see eqs. (B9)) was implemented for element TRI06_09. 
This stress field is represented by incomplete second-order 
polynomials; however, it identically satisfies the equations of 
equilibrium and does not possess spurious zero-energy modes. 

Four-node quadrilateral elements: QUA04_05, QTJA04J)7 y 
and QUA04_12. — Four-node elements have eight displace- 
ment degrees of freedom; thus five forces are necessary in the 
stress field approximation. The five-force field given in equa- 
tions (B8) was implemented for element QUAQ4_05 and re- 
sulted in a statically determinate four-node element. A complete 
linear polynomial was implemented for element QUA04_07, 
and a complete quadratic polynomial was used for element 
QUA04_12. The equations of equilibrium were identically 
satisfied for all elements and no spurious zero-energy modes 
were detected. 

Eight-node quadrUateralelements: QUA08J3, QUA08JL5, 
and QUA08J.8 . — Because eight-node elements have 16 dis- 
placement degrees of freedom, they require at least 13 indepen- 
dent forces to approximate the stress field. Quadratic 
polynomials do not contain a sufficient number of terms, so 
cubic polynomials must be used. Complete cubic polynomials 
(see eqs. (B5» were used for element QUA08„1 8, and reduced 
polynomials were used for element QUA08_15. These ele- 
ments were insensitive to rotation of the coordinate axes and 
did not possess spurious zero-energy modes. The stress field 
suggested by Nagabhushanam (Indian Institute of Science, 
Bangalore, India, personal communication, 1992; see 
eqs. (BIO)) was also implemented for element QUA08_13. A 
quadratic field that did not satisfy the equations of equilibrium 
resulted. It also possessed two spurious zero-energy modes. 

For all elements presented here, numerical integration was 
used to calculate the element matrices. In the case of the 
triangular elements, one-point integration was used for element 
TRI03_03, the three-point rule was used for elements with 
linear interpolation of the stress field, and the seven-point rule 


was used for elements with quadratic interpolation of the stress 
field. The locations of integration points were taken from 
reference 28. In the case of quadrilateral elements, standard 
Gauss integration was employed, with the 2 x 2 rule for 
elements with linear approximations of the stress field, the 
3x3 rule for elements with quadratic approximations, and the 
4x4 rule for elements with cubic approximations of the stress 
field. 

A patch test was performed for the present elements. Stress 
boundary conditions were prescribed for a finite element model 
of the test problem taken from reference 29, and all elements 
from the present library passed the patch test. 


Numerical Examples 

A number of example problems are presented in this section. 
Extensive numerical experiments were performed in order to 
establish the validity and accuracy of the Integrated Force 
Method, as well as to assess the relative performance of the 
present elements. The results obtained with present develop- 
ments are compared herein with corresponding analytical solu- 
tions. For some problems the responses obtained from the 
standard displacement method are also given in order to assess 
the potential advantages of the Integrated Force Method. The 
eight-node isoparametric element (ref. 27) was used in all 
displacement method calculations. 

Example 1: Bending of a Uniform Cantilever Beam 

Consider a cantilever beam of length L and uniform rectan- 
gular cross section dby H, as shown in figure 2(a). Assume that 
the beam is subjected to two distinct load cases: (1) a concen- 
trated force of intensity P, and (2) a uniformly distributed load 
of intensity q, and assume that the beam is made of a homoge- 
neous and isotropic material with a modulus of elasticity E and 
Poisson’s ratio v. By using two-dimensional finite element 
discretizations and assuming a state of plane stress, we can 
analyze the response of the beam. For this case the entire element 
library was implemented in order to establish the relative 
performance of the present elements. 

The influence of element shapes on the results was also 
investigated. Finite element discretizations obtained by using 
quadrilateral- and triangular-shaped elements are shown in 
figure 2, parts (b) and (c), respectively. The support conditions 
for modeling the beam with a clamped end assumed point a to 
be completely fixed and the horizontal displacements at points 
b and c to be suppressed. The circles in parts (b) and (c) of fig- 
ure 2 denote comer nodes, and the asterisks denote midside 
nodes. In discretizations using three-node triangular elements 
and four-node quadrilateral elements, midside nodes are not 
present. The dashed lines in figure 2(b) represent quadrilateral 
elements of distorted shapes. The concentrated force P was 
modeled by using nodal forces of intensities Pi and P 2 and 
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Figure 2.— A uniform cantilever beam showing geometric char- 
acteristics, loadings, and two-dimensional finite element 
models, (a) Subjected to concentrated force p and uniformly 
distributed load q. (b) Discretization using quadrilateral 
elements, (c) Discretization using triangular elements. 


assuming a parabolic distribution of the shear stress along the 
free end; the results were Pi = 0.5P, for elements with linear 
interpolation of geometry; and Pi = 0.1P and P 2 = 0.8P, for 
elements with quadratic interpolation of geometry . 

To analyze the beam for load case (1), convergence of the tip 
displacement was studied for discretizations by using various 
numbers of elements; the results are shown in figure 3 for 
triangular elements and in figure 4 for quadrilateral elements. 
The tip displacements were normalized with respect to the exact 
solution from the beam theory, including the average effect of 
shear stresses. Figure 3 shows that three-node triangular elements 
lead to very slow convergence of tip displacement, whereas six- 
node triangles provide accurate results with a relatively small 
number of elements and corresponding independent forces. 
The results shown in figure 3 also reveal that increasing the 
order of the stress approximation for three-node triangles does 
not improve the accuracy. 

Tip displacement convergence of quadrilateral elements was 
first studied with discretizations using elements of rectangular 
shape. The results, presented in figure 4(a) , show that all eight- 
node elements provided accurate results with a relatively small 
number of independent forces. Element QUA08J3, however, 
led to spurious zero-energy modes. Some additional degrees of 
freedom had to be suppressed in order to obtain a stable struc- 
ture. For the present analysis, vertical displacements at points 
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Figure 3.— Convergence study of tip displacement of 
cantilever beam using triangular elements. 


b and c were set to zero. These additional restraints resulted in 
convergence in the opposite direction for element QUA08_13. 
Figure 4(a) also shows the results for four-node quadrilateral 
elements. Element QUA04_05 provided a fast convergence, 
whereas elements QUA04_07 and QUA04_12, which use 
higher order stress approximations, produced stiff structures. 
Note, however, that both of these elements provide results 
within 0.5-percent error for a 24-element model. 

The effect that distortion of the element shapes has on the 
results was also studied. The distorted meshes were obtained by 
moving the comer nodes a distance of l m = 0.21^, as shown in 
figure 2(b), where L m is the length of the corresponding 
rectangular element. The results (see fig. 4(b)) show that four- 
node elements, especially element QUA04_05, are signifi- 
cantly less accurate. The eight-node elements, however, are 
almost insensitive to distortion in the model. 

Now let us consider stress distributions on the beam for load 
case (2). The intensity of the distributed load is taken as 
q = 12 kN/m; the length as L .= 12 m; the cross section 
dimensions as d == # = 1 .0 m; and Poisson 1 s ratio as v = 0.3. The 
values for normal stresses a and shear stresses r along the line 
y — -y s = -0.2887 m, which were obtained by using eight-node 
quadrilateral elements, are shown in figure 5 along with those 
calculated by using the beam theory. The locations shown in 
figure 5 coincide with the Gauss points for the 2 x 2 integration 
rale, which have been shown to be optimal sampling points for 
stress calculations in the displacement and hybrid methods 
(ref. 4). Figure 5 shows there is good agreement between the 
present results and the analytical solution and that results 
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Tip displacement, ^ v exact T >P displacement, \^ v exact 




Figure 4.— Convergence study of tip displacement of 
cantilever beam using quadrilateral elements, (a) Regular 
meshes, (b) Distorted meshes. 
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Figure 5,— Stress distribution along line y = “y g of cantilever 
beam subjected to uniformly distributed load, (a) Normal 
stress, (b) Shear stress. 


obtained with element QUA08_13 are in excellent agreement 
with the corresponding analytical solutions. However, caution 
must be exercised in using this element because of the presence 
of spurious zero-energy modes. 

Example 2: Pure Bending of a Circular Arch 

A circular arch of radius r 0 and rectangular cross section d by 
H, as shown in figure 6(a), is considered next. The arch is 
assumed to be clamped at 6- 0°, loaded with a concentrated 
moment of intensity Mat 0= 90° (where 0 is an angular 
coordinate), and made of a homogeneous and isotropic mate- 
rial with parameters E and v. This example is presented with the 
specific purpose of demonstrating the validity of the present 
elements in modeling domains with curved boundaries. A state 
of plane stress was assumed, and the arch was modeled by using 
two-dimensional finite element discretizations, as shown in 
figure 6(b). The boundary conditions for the clamped end were 
the same as those applied to the cantilever beam in the previous 
example. The circles and asterisks shown in figure 6(b) have the 
same meaning as in Example 1 . The nodes denoted by asterisks 
were not present when the arch was discretized with four-node 
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Figure 6-— A circular arch subjected to concentrated 
moment m . (a) One-dimensional model of a circular 
arch, (b) Two-dimensional finite element 
discretization. 

quadrilateral elements. In the finite element discretization, 
radii r a and denote the inner and outer contour, respectively, 
with r a = r 0 - 0.5 H, and r b = r Q + 0.5 H. The concentrated 
moment M was modeled by using the concentrated forces of 
intensities P\, and P3, which correspond to the exact stress 
distribution due to the concentrated moment. 

Let us first analyze the arch for r a = 10 m and = 12 m. Such 
an arch may be characterized as a thin arch (ref. 4). The remain- 
ing parameters are taken to be d = 1.0 m; E = 21xl0 7 kN/m 2 ; 
v = 0.3; and M == 600.0 kNm. The intensities of concentrated 
forces were calculated as P\ = 270.3 kN, P 2 = 56.3 kN, and 
P3 = 326.6 kN for discretizations with eight-node elements, and 
as Pi = P2 = 300 kN for discretizations with four-node 
elements. 

Results obtained with the present quadrilateral elements in a 
convergence study of the horizontal component u of the tip 
displacement, along with those obtained with an eight-node 


displacement isoparametric element, are shownin figure 7. The 
tip displacements were normalized with respect to the exact 
solution, which was calculated from the plane stress theory 
(ref. 23). The present elements, especially those with a quad- 
ratic interpolation of geometry, performed well. Note that the 
results for element QUA04__05 were obtained with an element 
local coordinate system such that the local x-axis is defined by 
the element centroid and the center of one of the element sides . 
The results obtained with the local axes parallel to the global 
axes are not shown because the responses exhibited unstable 
oscillations. This behavior, which is due to representing the 
stress field in terms of incomplete polynomials, reveals the high 
sensitivity of element QUA04_05 to the orientation of local 
coordinate systems. It also demonstrates the benefits of 
employing stress field representations composed of complete 
polynomials in the analysis of general two-dimensional 
problems. 

The displacements u along the line r~r Q were calculated next 
by using the discretization with six eight-node elements (see 
fig. 8). Exact displacements, calculated from the beam theory, 
are also shown in figure 8 for comparison. Again, the results are 
in good agreement. 

Stress distributions for a circular arch were also calculated. 
The results for normal stresses <r r and a t along the line 
r = r g = 10.423 m are shown in figure 9; they are compared with 
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Figure 7. — Convergence study of tip displacement u of a 
thin circular arch. 
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Figure 8.— Displacement distribution along line r= 2^ of thin 
circular arch. 




Figure 9. — Stress distribution along line r= 3^ of thin circular 
arch, (a) Radial stress, (b) Tangential stress. 




Figure 10. — Stress distribution along line r= Xg of a thick 
circular arch, (a) Radial stress, (b) Tangential stress. 


corresponding solutions from the displacement formulation 
and with the exact solution (ref. 23). Elements QUA08_15 and 
QUA08_18 performed slightly better than the displacement 
element. Figure 9 also shows that element QUA08_1 3 does not 
provide good stress predictions in this case. The stress field 
used for this element was obviously constructed to exactly 
model the beam bending. However, it does not satisfy Navier’s 
equations of equilibrium, and it does contain spurious zero- 
energy modes for some configurations. These characteristics 
make element QUA08_13 unsuitable for modeling general 
problems of two-dimensional elasticity. 

Stress distributions were also calculated for an arch with 
dimensions r a = 1.0 m, = 2.0 m, r g = 1.211 m, and M = 
300.0 kNm, with the remaining parameters being the same as 
before. Such an arch can be characterized as a thick arch. The 
concentrated forces used to model the moment M for this case 
were P x = 338.2 kN, P 2 = -73. 1 kN, and P 3 =265.1 kN. Normal 
stresses a r and a t were calculated by using the present eight-node 
elements; the results, together with corresponding analytical 


10 



solutions, are shown in figure 10. Note that elements QUA08_15 
and QUA08_18 provide more accurate results than the 
isoparametric displacement element, especially for radial stress 
<7 r , whereas element QUA08__13 exhibits the same problems as 
were observed for the thin arch. Also note that the reduced 
number of independent forces in element QUA08_15, as com- 
pared to element QUA08_1 8, did not lead to a loss of accuracy. 

Example 3: A Rectangular Plate Under Sinusoidal Load 

A rectangular plate with dimensions 2 L by 2a is shown in 
figure 1 1 along with its support conditions. The plate is sub- 
jected to a vertical load of intensity q(x) = q 0 sin(7u/2L). Since 
the geometry and the loading of the plate are symmetric with 
respect to the vertical axis, only half of the plate need be 
analyzed. The plate’s finite element discretization using quad- 
rilateral elements is also shown in figure 11. A state of plane 
stress was assumed for this analysis. The numerical values of 
the parameters were taken as E = 21.0xl0 7 kN/m 2 , v = 0.3, 
L = 10.0 m, a = 5.0 m; and q Q = 10.0 kN/m 2 . For comparison, 
the analytical solution of the problem was derived by following 
the procedure outlined in reference 23; this derivation is given 
in appendix C. 

The stress response was calculated for locations lying along 
the line AB by using both four- and eight-node quadrilateral 
elements, that is, with a mesh of 6 x 6 four-node elements and 
4 x 4 eight-node elements, respectively. The stresses were first 
calculated for locations corresponding to the optimal sampling 
points for the displacement and hybrid methods. These loca- 
tions correspond to Gauss integration points for lxl point 
integration when four-node elements are used, and 2x2 point 
integration when eight-node elements are used (ref. 24). Fig- 
ures 12 and 13 show the results for four- and eight-node 


q(x) 


v=Q 


romnrfc^ 




\ 

\ 

\ 

\ 


/ 

/ 



\ 

V 


/ 

/ 

/ 

/ 

L 




n x n elements 


v=0| 


Figure 1 1 .—A rectangular plate under sinusoidal loading. 


elements, respectively. These results agree well with the ana- 
lytical solutions given in equations (Cl), especially those for 
the eight-node elements. The present elements performed 
slightly better than the displacement formulation. 

The influence of the location of sampling points on the 
accuracy of stress predictions was also investigated. The 
stresses were calculated by using eight-node elements at loca- 
tions corresponding to 3 x 3 Gauss integration points. The 
results, shown in figure 14, indicate that the present elements 
are less sensitive to the location of sampling points than the 
corresponding displacement elements. It may also be con- 
cluded that the present elements provide better overall stress 
approximations within the element domains. 




Figure 12. — Stress distribution along line ab of a rectangu- 
lar plate as determined by using four-node elements at 
1 -point Gauss integration locations, (a) Normal stress a ■*. 
(b) Normal stress <r y (c) Shear stress 
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Figure 1 3. — Stress distribution along line ab of a rectangu- 
lar plate as determined by using eight-node elements at 
2-point Gauss integration locations, (a) Normal stress cr^ 
(b) Normal stress oy. (c) Shear stress 


Displacement calculations for the plate were also performed. 
The horizontal component u of displacement was calculated at 
element nodes lying along the line x=O.In figure 15 the results 
are compared with the exact solution given in equations (C2). 
A good agreement of the results is shown. 

Example 4: A Rectangular Plate With a Circular Hole 
Under Uniform Tension 

A rectangular plate of dimensions 2 a by 2b with a circular 
hole of radius r, as shown in figure 16(a), was analyzed. The 
plate was assumed to be subjected to a uniform tension of 
intensity q along the Ox-axis and made of a homogeneous and 
isotropic material. Because of the double symmetry of the 


0 2 4 6 8 10 

(a) Distance along x-axis, m 



Figure 14. — Stress distribution along line ab of a rectangu- 
lar plate as determined by using eight-node elements at 
3-point Gauss integration locations, (a) Normal stress 
(b) Normal stress <r T (c) Shear stress 


geometry and the loading of the plate, only a quarter of the plate 
(i.e., that bound by the arc AB and lines BC, CD , DE, and EA) 
was discretized with finite elements. Finite element 
discretizations using quadrilateral and triangular elements are 
shown in figure 16, parts (b) and (c), respectively. The symme- 
tryboundary conditions u - 0 along the line AE, and v = 0 along 
the line BC were applied. A state of plane stress was assumed 
in the analysis. 

The stress concentration factor at point A of the plate was 
calculated for a = 48 cm, £ = 24 cm, and r = 6 cm by using the 
entire element library. The results, given in table I, are com- 
pared with the adjusted stress concentration factor calculated 
by using the expression from reference 30 (i.e., k = 3.2126 for 
the given dimensions of the plate). A very good agreement for 
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0.8x10-6 



Distance along y-axis, m 

Figure 1 5.— Horizontal displacement along line x= 0 of 
rectangular plate. 


the stress concentration factor was achieved with the elements 
from the present library, especially with those with quadratic 
interpolation of geometry, as depicted in table L This example 
further demonstrates the accuracy of the Integrated Force 
Method in stress calculations, particularly at locations that do 
not coincide with optimal sampling points. 

Discussion 

The validity and accuracy of the finite element library 
presented in this paper were demonstrated through numerical 
examples. Both one- and two-dimensional problems of elastic- 
ity were analyzed, and from the numerical results previously 
presented, a comparison can be made of the element perform- 
ances. A careful examination of numerical results reveals that 
elements with quadratic interpolation of geometry performed 
better than those with linear interpolation of geometry . This is 
only partly due to the higher order of approximation of the 
stress fields, and thus, larger number of independent forces. 
Such a conclusion is supported by the fact that an increase in the 
number of independent forces in three-node triangular elements 
does not result in improved accuracy, and that higher order 
approximations of stress fields for four-node quadrilateral 
elements may result in overly rigid models. We may conclude 
that the interpolations of stress and displacement fields cannot 
be chosen arbitrarily, but must be compatible with stress-strain 
law. Moreover, elements QUA08_1 8 and QUA08_1 5 exhibited 
better overall performance than element QU A08_l 3 . The stress 
fields used for these two elements are represented by complete 





Figure 16. — Finite element discretizations of rectangular 
plate with circular hole in uniform tension, (a) Two- 
dimensional plate, (b) Discretization with quadrilateral 
elements, (c) Discretization with triangular elements. 


TABLE I. — STRESS CONCENTRATION FACTOR AT 
LOCATION A FOR THE PLATE WITH A HOLE 


Element 

Number of 

Number of force 

Stress concentration 

type 

elements 

unknowns 

factor 

TRI03_03 

37 

111 

2.54649 

TRI03_05 


185 

2.54649 

TRI03_07 


259 

2.54649 

TRI06J39 

37 

333 

3.25875 

TRI06_1 1 


407 

3.15126 

TRI06_12 


444 

3.10322 

QUA04 05 

30 

150 

2.74766 

QUA04_07 


210 

2.89305 

QUA04_12 


360 

2.96215 

QUA08.13 

30 

390 

3.00314 

QUA08_15 


450 

3.25266 

QUAG8_18 


540 

3.26137 
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third-order polynomials, which produce elements that are in- 
variant with respect to transformation of the local coordinate 
systems and that are free of spurious zero-energy modes. These 
stress fields also identically satisfy the equations of equilib- 
rium. The stress field used for element QUA08„1 3 , on the other 
hand, was represented by second-order polynomials, which 
yielded two spurious zero-energy modes for the rectangular 
configuration of the element. In addition, it did not satisfy the 
equations of equilibrium, and thus led to erroneous stress 
calculations in Example 2. 

For some problems, good results were obtained with element 
QUA04J35, which has a very small number of unknown 
independent forces. The tip displacement convergence study in 
Example 2 revealed, however, that this element is sensitive to 
the orientation of coordinate axes and may not always be 
reliable in the analysis of domains with arbitrary geometric 
configurations. Moreover, element QUA04_05 provides accu- 
rate stresses only in the centroid, which may not always suffice. 
Elements with quadratic interpolation of geometry produce 
more accurate stress predictions at the optimal sampling points 
and at arbitrary locations within the element. Similar conclu- 
sions can be drawn for triangular elements. 

The Integrated Force Method was also compared with the 
assumed displacement based finite element method. The re- 
sults presented here reveal that the Integrated Force Method is 
better overall for stress calculations and provides displacement 
predictions of comparable accuracy. 

Concluding Remarks 

A finite element library was developed to analyze two- 
dimensional structural mechanics problems by the Integrated 
Force Method. Triangular- and quadrilateral-shaped elements 
capable of modeling domains with arbitrary geometric con- 
figurations were presented. The displacement and stress fields 
were independently approximated. Displacement interpolation 
was performed as in the standard displacement method, and a 
procedure was developed to derive the stress interpolation 


functions in terms of complete polynomials of the required 
order. An Airy stress function was written as a complete 
polynomial of order/? that contains p + 1 independent constants. 
The definitions of stress components in terms of stress func- 
tions were used next to derive the expressions for stresses. 
Elimination of the dependent constants from the expressions 
for stresses yielded stress fields expressed in terms of complete 
polynomials of order p - 2. Stress fields thus defined identically 
satisfied the equations of equilibrium. The resulting element 
matrices had the correct rank and were insensitive to the 
transformation of local coordinate systems. 

The present elements were applied to solve a variety of 
problems in two-dimensional elasticity. Comparisons were 
made with corresponding analytical solutions, and there was 
good agreement of the results. A series of numerical tests were 
performed in order to assess the relative performances of the 
present elements. These studies showed that elements with 
quadratic interpolations of geometry and displacements pro- 
vide reliable predictions for all problems. The four-node quad- 
rilateral elements performed well for some problems. Element 
QUA04__05 was shown to provide good results for a small 
number of unknown quantities, but it was sensitive to the 
orientation of the local coordinate system. This condition 
restricts its range of application. 

The Integrated Force Method was also compared with the 
standard displacement method. The results presented here 
reveal that overall the Integrated Force Method performs better 
in stress calculations and exhibits an accuracy in displacement 
predictions comparable to the displacement method. These 
results confirm that the Integrated Force Method can be used 
successfully and efficiently in structural analysis and provide 
justification for efforts to incorporate the force method of 
analysis into general purpose finite element programs. 


Lewis Research Center 
Cleveland, Ohio 
July 17, 1995 
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Appendix A 
Symbols 


A c 

complementary energy 

m 

number of system independent forces 

A P 

potential energy 

m e 

number of element independent forces 

a , b 

plate dimensions 

[N] 

matrix of displacement interpolation functions 

[BJ 

element equilibrium matrix 

Ns 

number of prescribed displacements 

[BJ 

portion of system equilibrium matrix correspond- 
ing to modes with prescribed displacement bound- 

N t 

total number of displacement degrees of freedom 


ary conditions 

n 

number of system equilibrium equations 

[BJ 

portion of system equilibrium matrix correspond- 
ing to modes where external loads are prescribed 

n e 

number of element displacement degrees of 
freedom 

[C] 

compatibility matrix 

n r 

rank of element equilibrium matrix 

9 

arbitrary constants 


intensity of concentrated force 

[D] 

compliance matrix of material 

{P} ' 

system equivalent load vector 

d 

cross-sectional dimension 

{P*} 

total of right side of system of equations 

E 

modulus of elasticity 

{P e } 

vector of element equivalent nodal loads 

{¥} 

system vector of independent forces 

P 

order of polynomial 

{F e } 

vector of element independent forces 

Mo 

intensity of distributed load 

Fi 

generalized force coefficients 

{R} 

vector of support reactions 

[G] 

system flexibility matrix 

r 

number of compatibility conditions 

[G e ] 

element flexibility matrix 

r a> r b> r o 

radial coordinates 

H 

cross-sectional dimension 

[S] 

system matrix 

[J] 

n x m deformation matrix 

{U} 

vector of system unknown nodal displacements 


length 

{U e } 

vector of element nodal displacements 

[L] 

matrix of differential operators that defines strain 
displacement relationship 

{«} 

displacement vector 

l 

number of element rigid body modes 

u,v 

displacement components 

i 

distance comer node is moved 

V 

volume 

M 

intensity of moment 

x,y 

Cartesian coordinates 
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m 

stress interpolation matrix 


components of strain vector 


coordinate of Gauss point 

e 

angular coordinate 

m 

[L] [N] 

V 

Poisson’s ratio 

ip,} 

vector of element deformations 

{ 0 } 

stress vector 

{P o) 

vector of initial deformations 

G x Gy 

components of stress vector 

y*y 

component of strain vector 

% 

component of stress vector 

{6 0 } 

effective deformation vector 

0 

linear functions of constants 

{£} 

strain vector 

0 

Airy stress function 
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Appendix B 

Expressions for Stress Fields 


(a) Full Polynomials Derived from the Stress Function 


- constant terms: 

= f/ 0) 

(Bla) 


4 0) - f 2 (0) 

(Bib) 


^ = # 

(Blc) 

- linear terms: 

4 1} = Fi (1) x + F 2 (1) y 

(B2a) 


4 J) = F^x+F^y 

(B2b) 


r£> = -F^x-F^y 

(B2c) 

- quadratic terms: 




= - iF 5 ( V + F<V + F 2 2) xy 

(B3a) 

4 2) 

= Fi 2) x 2 -iFi 2) y 2 + Fi 2) xy 

(B3b) 

4 ? 

= -^i 2) * 2 -iF 2 (2) y 2 + F 5 (2) xy 

(B3c) 

- cubic terms: 



= 

-JjtfV + FpV - Fj 3) x 2 y + F^xy 2 

(B4a) 

a y — 

fJ 3) x 3 - >y 3 + Fi 3) x 2 y - i^ 3) xy 2 

(B4b) 

r ( 3 ) - 
>xy 

-|fJ 3) - |f| 3) + F 5 (3) x 2 y + Fj 3) xy 2 

(B4c) 

- full cubic polynomial 



Ox = Fi + F 4 x + Fsy - |Fi 2 * 2 + Fgy 2 + Fgzy - 


\Fnx 3 + F l3 y 3 - F 18 x 2 y + F w xy 2 

(B5a) 


cr y = F 2 + F 6 x + F r y - F w x 2 + ^F i2 y 2 + F n xy + 

Fisx 3 - ^Fisy 3 + Fi6* 2 y - F u xy 2 (B5b) 

Txy = F 3 - F 7 x - F 4 y - |F u a: 2 - ^F 9 y 2 + Fi 2 xy - 

|fi 6 x 3 - |Fi 4 y 3 + Fnx 2 y + F ls xy 2 (B5c) 
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(b) Full Polynomials Reduced Using the Condition V 2 (crj. + a y ) = 0 

- quadratic terms: 

4 2r) = F< 2r) (y 2 - l* 2 ) + if r) xy - ijfV (B6a) 

4 2r) = -fif’V + rt } (* 2 - ^y 2 ) + Fp r) *y (B6b) 

rg r ) = Ff r, *y-|F 2 (2r) y 2 + F3 (2r) xy-|F 4 (2r) x 2 (B6c) 

- cubic terms: 

•fJ = if* (V - J* 2 y) + if* (*»* - £* 3 ) - iff’* 3 - ^ 3r) * ! V (B7a) 

<f> = + (B7b) 

= |j^ 3r V + fi 3r, (i» 3 »-|» 3 )+|4’ r) * 2 »+d ar) (5*» , -5* 3 ) (B7c) 

(c) Stress Field Used for the Element QUA04_05 

°x = Fi + F 4 y (B8a) 

<T y = F2 + F5X (B8b) 

r*,, = F 3 (B8c) 

(d) Stress Field Used for the Element TRI06.09 

<t x = Fj + F 2 y + F 6 x - 2F 8 xy (B9a) 

<r y — F 3 + F 4 x + Fry - 2Fgxy (B9b) 

T xy — F 5 - F 6 y - Frx + F 8 y 2 + F 9 x 2 (B9c) 

(e) Stress Field Used for the Element QUA08.13 

o x — Fi + F 2 y 4- F 6 x + F 8 y 2 + Flo* 2 - 2F 13 xy (BlOa) 

<7 V = F 3 + F 4 x + F 7 y + F 9 x 2 + Fuy 2 - 2Fi 2 xy (BlOb) 

Txy = F5-F6y-Frx-2(F 10 + Fu)xyH-Fi 2 x 2 + Fi 3 y 2 (BlOc) 
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Appendix C 

Analytical Solution for the Rectangular Plate Under Sinusoidal Load 


The analytical solution of the problem presented in Exam- 
ple 3 is given here. It can be derived by using the Fourier series 
outlined in reference 23. The stress components a x , a y , and % ^ 
are given as 


2 (2 

a x ~a Cj cosh ay + C 2 sinh oy+C 3 — sinh ay + y cosh ay 


where C\ = /^oCsinh aa + aa cosh aa)/D\, C 2 ~ -po(cosh oca + 
aa sinh cca)/D 2 , C 3 = p$a cosh aalD 2 ; and C 4 = -p 0 a sinh oaf 
Pi\ with D\ = a 2 (2aa + sinh 2aa ), D 2 = 2oa ^ sinh 2a<z), 

a= tt/2L, and x and y are coordinates of the point inside the plate 
as defined in figure 1 1 . 

The displacement w is calculated by integrating the strain- 
displacement relations as follows: 


f 2 \ ^ 

+ C 4 — cosh oy + y sinh ay sin — (Cla) 
\a ) 2 L 


cc 

u = - — |Cj (1 + v) cosh ay + C 2 (1 + v) sinh ay + C 3 


CT y = ~a 2 (Cj cosh qy + C ' 2 sinh qy + C$y cosh qy 


+C 4 jcosh qy) sin — (Clb) 


-sinh qy + (l + v)ycoshqy 


( 2 'N yjX 

+ C 4 I— cosh oy + (1 + v)y sinh ay cos— (C2) 
\a J 2L 


2 /I 

= a; sinh ay + C 2 cosh ay + C 3 [ — cosh ay + y sinh ay 

Z 1 \ 

+ C 4 —sinh ay + y cosh ay cos — (Clc) 
v a J 2 L 
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